function b = adsecurity_hhw(Xit,Cg,MU,Y,CAPH,T,rho,Nstates)
%This function gives the ad security for the future period given history.
%The output is b=[b(lambda=1|Lambda_t); b(lambda=2|Lambda_t),...,b(lambda=n|Lambda_t)]

%In the first step we get the current period, spot price and guarantted
%price from history Lambda_t
    t0=length(Xit)+1;
    Y1 = repmat(reshape(Y,1,1,T),Nstates+1,1);
    C0=Y1-MU;
    C0(Nstates+1,:,:)=0;

    c0=C0(:,:,t0);
    tildec=findtildec(Xit,Cg,Y,MU,Nstates,T);
    barct=Cg(:,t0);
    n=Nstates+1;
    tildec=max([barct repmat(tildec,n)],[],2);

    if t0==T
        %if current period is final period, 
        % The bond is guaranteed consumption minus endowed consumption
        b=tildec-c0;
        b(b<0)=0;
    else
        %If not, we first get all future spot consumption.
        tf=T-t0;
        barcf=Cg(:,(t0+1):T);
        % future endowment
        c0f = C0(:,:,t0+1:T);

    
        %Then we calcualte the aggregate probability agent would not lapse
        %(same as denominator of equation 1 in the paper)
        %First, we get the state where the agent would stay
    
        Nstates = size(CAPH,1);
        barcf=reshape(barcf,[Nstates 1 tf]);
        stay  = repmat(permute(barcf,[2,1,3]),Nstates,1)<repmat(tildec,1,Nstates,tf);
    
        % death is a no stay state
        stay(:,Nstates,:) = 0;
        % probability of getting to each state in a 
        % given period -- transit from periods where i stay 
        capH=CAPH(:,:,t0);
        pt = zeros(Nstates,Nstates,tf);
        pt(:,:,1) = capH;
        
        E = zeros(Nstates,Nstates,tf);
    
        E(:,:,1) = rho*pt(:,:,1).*stay(:,:,1).*repmat(permute(c0f(:,:,1),[2,1,3]),Nstates,1,1);

        D = zeros(Nstates,Nstates,tf);
        D(:,:,1)= rho*pt(:,:,1).*stay(:,:,1);
        for t=2:tf
            capH=CAPH(:,:,t0+t-1);
            pt(:,:,t) = pt(:,:,t-1).*stay(:,:,t-1)*capH;
            E(:,:,t) = rho^t*pt(:,:,t).*stay(:,:,t).*repmat(permute(c0f(:,:,t),[2,1,3]),Nstates,1,1);
            D(:,:,t)= rho^t*pt(:,:,t).*stay(:,:,t);
        end
        % Add them up

        % numerator
        E = sum(sum(E,2),3);
        num = c0+E;
        %denumerator
        den = 1 + sum(sum(D,2),3);
        %use equation 3 in ra note to calculate bond.

        b=den.*tildec-num;
        b(Nstates)=0;
    end
end